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Abstract 

This work describes a new ID hybrid approach for modeling atmospheric pressure discharges fea¬ 
turing complex chemistry. In this approach electrons are described fully kinetically using Particle- 
In-Cell/Monte-Carlo (PIC/MCC) scheme, whereas the heavy species are modeled within a fluid 
description. Validity of the popular drift-diffusion approximation is verified against a ’’full” fluid 
model accounting for the ion inertia and a fully kinetic PIC/MCC code for ions as well as electrons. 
The fluid models require knowledge of the momentum exchange frequency and dependence of the 
ion mobilities on the electric field when the ions are in equilibrium with the latter. To this end 
an auxiliary Monte-Carlo scheme is constructed. It is demonstrated that the drift-diffusion ap¬ 
proximation can overestimate ion transport in simulations of RF-driven discharges with heavy ion 
species operated in the 7 mode at the atmospheric pressure or in all discharge simulations for lower 
pressures. This can lead to exaggerated plasma densities and incorrect profiles provided by the 
drift-diffusion models. Therefore, the hybrid code version featuring the full ion fluid model should 
be favored against the more popular drfit-diffusion model, noting that the suggested numerical 
scheme for the former model implies only a small additional computational cost. 
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I. INTRODUCTION 


The discharges operated under atmospheric pressure are easier to operate compared to the 
low-pressure discharges because the expensive vacuum equipment is no longer needed. How¬ 
ever, due to their smaller size compared to the low-pressure counterparts, the experimental 
diagnostics proves to be more difficult. Hence, numerical modeling plays an important role 
in facilitating the understanding of processes taking place in atmospheric pressure plasma 
discharges. 

If the scale of spatial inhomogeities in a plasma discharge is much greater than the elec¬ 
tron energy relaxation length A^- ~ A( 2 me/M)^/^ S> A (with A the momentum exchange 
mean free path and the characteristic time scale of phenomena in interest is larger than 
the corresponding time scale of the energy relaxation = l/((5z/)), the spatial derivatives 
in the Boltzmann’s equation governing the electron energy distribution function (EEDF) 
can be neglected and the distribution function will only implicitly depend on the spatial 
inhomogeneities through the spatial dependence of the electric field and the plasma param¬ 
eters. In this case one can use local description for plasma behavior. However, since in most 
atmospheric plasmas Coulomb collisions ” maxwellizing” the EEDF are negligible compared 
to the electron-neutral collisions (usually leading to pronouncedly non-maxwellian EEDFs) 
one still has to determine the EEDF and calculate the corresponding transport coefficients 
to be used in a fluid model for electrons. Typically, to this end a Od Boltzmann solver nn 
is used. 

Whereas this popular approach is suitable for studying most atmospheric pressure dis¬ 
charges, there are situations where it turns out to be inadequate even in such highly col- 
lisional plasmas (e.g., ISHHI). In [5] it was shown that one of the indispensable ingredients 
in the mechanism of striation generation in an atmospheric dielectric barrier discharge is 
the non-local energy transport of the electrons. In [B] it was demonstrated for a glow dis¬ 
charge in helium in the 7 mode that in contrast to the expected two-temperature EEDF it 
exhibits three different energy groups of electrons. In addition to the non-Maxwellian ener¬ 
getic tail on the distribution function of the electrons and the mid-energy group of electrons 
common to the glow discharges operated in the 7 mode there was observed a low energy 
group of electrons with an unusually low temperature contributing most to the electron den¬ 
sity. Such a three-temperature EEDF is impossible to model with a fluid description which 
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usually assumes a maxwellian distribution function for electrons. The EEDF becomes even 
more complicated and farther from the maxwellian shape when additional species, especially 
molecular gases, are added to the discharge. Furthermore, recently it has been demonstrated 
in [7] that due to the large spatial nonuniformities of the electric field atmospheric pressure 
plasma discharges can exhibit distinctly nonlocal energy transport even in the hi (ohmic) 
mode. These effects obviously cannot be captured by the local description. Yet another 
example of nonlocal electron energy transport in the highly collisional plasmas is studied 
in [H], where an efficient gas impurity detector based on a ’’short” DC discharge, which 
predominantly consists of the negative glow region, is proposed. Since the transverse size of 
the discharge there is chosen to be smaller than the energy relaxation length, the energetic 
electrons born in Penning ionization reaction can diffuse to the discharge wall where the de¬ 
tector is placed, retaining a large portion of its energy and thus producing pronounced peaks 
in the EEDF. The energy of such electrons is T = T* — T/ with T* the excitation energy of 
the working gas (chosen to be helium for the high potential energy of its metastable levels) 
metastable atoms or molecules and Sj the ionization energy for different background gas 
components. Since Sj is distinct for different gas impurities, it is possible to identify the gas 
chemical composition by tracking peaks on the EEDF. If the electron energy transport had 
been local, the EEDF of the electrons leaving the discharge would have had the maxwellian 
shape with no peaks on it. 

Of particular importance is an appropriate description of energetic electrons. Although 
the number of such electrons is too small to affect the plasma bulk dynamics, these electrons 
are essential for the plasma sustainment as they signihcantly contribute to the particle 
balance through the ionization processes, which have threshold at relatively high energies. 
On the other hand, the lower energy electrons determine the bulk plasma dynamics, which in 
turn governs generation of the electric field accelerating the high energy electrons. Therefore, 
it is important to correctly reproduce the correct energy distribution function of the electrons 
describing all its groups. 

The details of ion energy distribution function in the atmospheric pressure plasma dis¬ 
charges are, in contrast, rarely of signihcance as they have much smaller collisional mean 
free path than electrons and thus have the Maxwellian form in all cases of practical impor¬ 
tance. This has motivated us to use a hybrid numerical scheme where ions are described 
using the fluid approximation and the electrons are followed kinetically. Furthermore, most 
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of the experimental data concerning the ion-neutral collisions are provided in the form of the 
reaction rates (it is much easier to hnd the differential cross-sections for the electron-atom or 
electron-molecule collisions) and not the energy-resolved collision cross-sections, which fur¬ 
ther renders the kinetic description of the ions an overkill. The fluid description employed 
for ions also simplihes bookkeeping of the potentially complex chemistry taking place be¬ 
tween the neutral and charged species and the kinetic description of the electrons ensures an 
accurate energy and particle balance in the numerical description of the discharge. An addi¬ 
tional advantage of using a fluid description for the heavy particle species is that it requires 
much less computer memory compared to the kinetic description. Although this argument 
does not play a signihcant role in ID simulations, it can become important as long as one 
wants to study multidimensional models with the kinetic part of the algorithm parallelized 
on graphics processing units (GPUs), as the latter usually have dedicated limited memory 
which cannot be expanded. 

In the present work we propose a hybrid approach, where electrons are followed utilizing 
the fully kinetic particle-in-cell technique and ions are described using a fluid model. Before 
one can trust the numerical results of a particular model, such a model has to be validated 
and verihed by comparing its results with either an analytical treatment or results of a well 
established model. Should such a comparision demonstrate any discrepancies, they should 
be amenable to an explanation based on the model differences. That is the main goal of the 
present work. 


II. DETAILS OF THE HYBRID CODE 

In the present work we study plasma dynamics across the discharge channel in a RF- 
driven plasma (CCP) planar discharge at atmospheric pressure, such a micro-jet discharge. 
Correspondigly, we limit ourselves to ID spatial dimension with coordinate being the 
distance from the powered electrode. The driven (grounded) electrode is located at ^ = 
Q [z = L). Following the logic outlined in the introduction section, we adopt a kinetic 
description for the electron and a fluid description for the ion and neutral species. The 
resulting approach enables an efficient implementation of numerical model for discharges 
featuring a complicated chemistry (e.g., m) 
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A. Kinetic Description of Electrons 


The electron dynamics is traced numerically using the PIC technique mm, where the 
electron particle distribution function is discretized on a moving Lagrangian grid with a 
number of markers following the characteristics of the Boltzmann’s equation, and the electric 
field is discretized on a stationary Eulerian grid. The markers used for the representation 
of the particle distribution function can be also regarded as ’’superparticles” representing a 
large number of physical particles, which are close to each other in phase space, so that 


fe{z,v) ='^WjW{Zj{t) -z)6{Vj{t) -v), (1) 

j 

where W is the superparticle shape function and Wj is weight of the jth particle. In case 
of a field grid uniform in z and homogeneous initial plasma density (which we take to be 
the case in the present work), it is convenient to take the same weight for all superparticles, 
Wj = Ueo/NgQ with Ugo the initial electron density and Ngo the initial number of superparticles 
per held grid cell. This results in the following number of real particles per superparticle, 

_ UeoAV 

^^p/sp — 

tVsO 


with AV the grid cell volume, AV = AzS, Az being the held grid cell size and S the 
electrode area, respectively. 

The PIC/MCC numerical scheme describing the electrons in our implementation con¬ 
sists of four parts: particle pusher, particle removal, charge density assignment, and the 
Monte-Carlo collisions module. In the particle pusher particles are moved in phase space by 
integrating the Newton’s equations of particle motion under inhuence of the electric held. 
This is done using the explicit leapfrog scheme with the electric held interpolated from the 
held grid nodes to the particle position with the interpolation function S' being the same 
as the particle shape function S in Eq.([^, which we chose to be of the Cloud-In-Cell (CIC) 
variety (e.g., M)- The superparticle weight remains constant during the particle pusher part 
of the PIC/MCC algorithm due to Liouville’s theorem. The pusher algorithm also checks if a 
particle gets rehected or generates new particles at the borders of the computational domain 
due to the corresponding surface processes, such as generation of the secondary electrons, 
by using corresponding probability. After the particle pusher is completed, it is checked if 
a particle leaves the computational domain. Should it be the case, the particle is removed. 
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Further, one calculates charge density of the electron component by extrapolating the su¬ 
perparticle charge to the surrounding held grid nodes. Finally, a Monte-Carlo technique is 
used to implement the collisions between the electron and the heavy species (charged and 
neutral). To this end, we exploit a variety of the ’’null-collision” method described in [12]. 
Rather than taking N^Pmax number of electrons for executing a collision algorithm with Ne 
the total number of electrons and Pmax = 1 — exp(—i/maxAt) the maximum ’’null-collision” 
frequency, one checks if i? < Pmax with R the pseudorandom number uniformly distributed 
in [0,1] for every particle. In electron collisions with the background neutral gas atoms or 
molecules the latter are treated as having a uniform density and any change of its density 
and temperature in the course of the discharge evolution is assumed to be negligible. In 
contrast, when treating collisions of electrons with the other species such as metastables and 
ions, the density prohles of the latter species are taken into account. When only a reaction 
rate is available, we deduce energy dependence from the reaction rate electron tempera¬ 
ture dependence based on the Maxwellian ansatz and then use the resulting energy-resolved 
cross-section in the simulations similar to [12]. A proper choice of pseudorandom number 
generator is also important as the number of collisions in simulating such highly collisional 
plasmas is large and a good pseudorandom generator should have an appropriately long 
period. For our hybrid code we used the 128 bit Xoorshift pseudorandom number generator 
suggested in [1^, which has period of 2^^® — 1. To accelerate the computations, the kinetic 
PIC/MCC part of the present hybrid code was implemented on a graphics processing unit 
(GPU) analagous to the GPU PIG/MGG code used in the benchmarking study [TB] . 

The cell size Az was limited by the need to resolve the Debye length in order to avoid 
excessive numerical heating (in our simulations we took Az = 0 . 2 Ai 5 e,max, where Xoe^max is 
Debye length calculated using the maximum expected electron density rie^max) and the time 
step was limited by the need to resolve the elastic collisions (we took At = with 

^max the maximum value of the elastic collisions over the computational energy domain). 

B. Fluid Description of Heavy Particle Species, Full Treatment 

In contrast to the kinetic treatment of the electron component, the heavy species transport 
is followed using a fluid model. In such a model the heavy particle density is governed by 
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the particle continuity equation, 


duhs 

dt 


+ V ■ F/is — Ghs I 


(3) 


where hs is the heavy particle species subscript, = UhsVLhs is the heavy particle flux with 
U/is the average ordered velocity of this species, and Ghs the net density change rate due to 
the reactions where such particles are either born or destroyed. In general case (see, e.g.. 
Cl) the heavy particle flux is governed by equation which can be obtained by multiplying 
the Boltzmann’s equation by particle velocity and integrating it over velocity space, which 
yields 

Mhs + V ■ ) = ZhseUhsE. - V{nhsThs) - ^hs,NMhsThs (4) 


with E the electric held, Zhs the species charge number, Mhs its mass, Ths the species 
temperature, Vhs,N the momentum exchange frequency of the species with the neutral back¬ 
ground gas iV, and the expression in the brackets of the second term on the left hand side is 
a dyadic tensor. In the atmospheric pressure plasma discharges the left hand side of Eq. (|^ 
is normally neglected and the ion fluxes are calculated from the right hand side of this 
equation (drift-diffusion approximation, see the next section). However, it is pointed out in 
nn that the right hand side of Eq. (|^ vanishes only after the equilibrium drift velocity is 
achieved and it takes several nanoseconds for a heavy ion to accomplish it. A large ion flux 
gradient, which can arise in a highly collisional discharge, for example, due to the intense 
ionizing electron avalanche in the 7 mode, can cause the second term on the l.h.s. of Eq. (|^ 
to be of significance as well. Since the drift-diffusion approximation is also frequently used 
for lower pressure discharges (up to 100 mTorr), the explicit time derivative in Eq. (|^ can 
clearly become comparable with the last term in this parameter range. All this suggests 
that the left hand side of the latter equation can be substantial and it is interesting to verify 
how close the results of the full model utilizing Eq. (|^ are to the results of the popular 
drift-diffusion approximation under different conditions. 

Finally, we note that the energy transport equation is usually omitted for the heavy 
species since they come into the thermodynamical equilibrium with the background neutral 
gas really fast due to the efficient collisional energy exchange in contrast to the electrons. 
Under this assumption the temperature is usually approximated as Tjv for the neutral species 


and by Eq. (10) for the ion species. 
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A convenient numerical scheme for solving Eqs. Q and ([^ employs the leapfrog time 
intergration approach, where discretization of the density n is performed at integral time 
levels, and the flux F is taken at time levels half a time step apart. In the particular scheme 
that we have used in this study one hrst calculates the particle flux at the next time level, 
so that in the one dimensional case the corresponding equations read 


rn+V2 A ^ 


_ 1 _ Af 
2 A 2 


/ 2 Y^n —1/2 

^j + 1 ^j + 1 


-i + l 


_ pn— 1/2 / 


^Ai 


■pTl-l- 1/2 -pTl—1/2 

^ j-1 ^j-1 

”i-i 


+ 


M 


At 

2AzM 


n 


n n~'n _ 'j^n \ 

j+i-Lj+i 


(5) 


for j = 1..N — 2, which can be easily found by integrating Eq. Q from Zj- 1/2 = {j — 
l/2)Az to Zj-^-i /2 = (j + 1 / 2 )A 2 :. The equations at the boundary nodes are obtained by the 
same integration procedure using the corresponding intervals [ 2 ^ 0 , 2 : 1 / 2 ] and [ 2 ;Ar_ 3 / 2 , 2 ;Ar_i], 
respectively. The procedure being straightforward, we omit the results. Note that the 
advection term on the l.h.s. of these equations is discretized semi-implicitly, which helps 
to keep the equations linear and tridiagonal. Then, they can be easily solved, for example, 
by using the Thomas’ method. Generally speaking, the momentum exchange frequency 
depends on velocity u = T/n, which is spatially nonuniform. To avoid the need for an 
iterative solver we use m” = which reduces the order of the numerical scheme, 

but we hnd in the benchmark comparisons with PIC simulations the resulting accuracy to 
be still sufficient due to the small time-step caused by the need to resolve the collisions in 
the PIC/MCC part of the algorithm describing electrons. A higher oder alternative would 
be to solve Eqs. (|^ and (|^ iteratively, which is more computationally expensive. Although 
the diffusion term proved to be small in all the cases we considered, we still retained it, 
albeit using a simplihed expression from the drift-diffusion approximation given in Eqs. @ 
and (10). To determine the velocity-dependent momentum exchange frequency u in Eq. ([^ 
we have constructed an auxiliary Monte-Carlo code, which calculates it based on the energy 
resolved collisional cross-sections. A detailed desciption of the latter code is given in the 
next subsection. 

After the particle flux at the next half-integer time level is found, one can obtain the 
value of the density at the next integer time level from the particle continuity equation. 

At 


= n" 


2Az 


('pn+ 1/2 _ pj+l/2j ^ A/G'^+1/2^ 


( 6 ) 


for j = 1..N — 2, which is derived by integrating the particle continuity equation over 
[ 2 ;/_i/ 2 , ; 2 /+i/ 2 ]- Equations for the bondary points j = 0 and j = N — 1 are obtained similarly 










by integrating over the corresponding intervals [zq, ^ 1 / 2 ] and [^Ar- 3 / 2 , and are not shown 

here. For the sake of simplicity in our numerical implementation we substitute 
with G”, which degrades accuracy of the scheme. However, by comparison with the PIC 
simulations we conclude that the resulting accuracy is still sufficient not to cause a significant 
deviation from the fully kinetic calculations. 


C. Fluid Description of Heavy Particle Species, Drift-Diffusion Approximation 


In the frequently used drift-diffusion approximation the particle fluxes are assumed to 
have reached a quasi-stationary value, so that the left hand side of Eq.(|^ is neglected. In 
this case the particle flux is approximated as 

T/is ZfigCTbhg^iigEi Dfig^Tifig (7) 


with 




the heavy particle mobility (note that it equals zero for neutral species), and D^g the heavy 
particle diffusion rate. 

The diffusion rate accounting for the impact of electric field on the ion diffusivity is 
obtained by using the generalized Einstein’s relation [19] , 

hiEj 


D, = 


Qi 


(9) 


with 

T. = Tn + ^ \ ^ ^ mME f- 10 

orrii + druN 

Using the expression for particle flux in the drift-diffusion approximation given in Eq. ([^ 
one can obtain particle density from Eq. ([^. This equation can be solved either explicitly (by 
taking the particle flux values not at the time level n-|-l/2, but n), in which case the boundary 
condition is not needed, or, for example, semi-implicitly (if the original form of Eq. (|^ is 
retained. In the latter case the boundary conditions come from the kinetic limitation on the 
particle flux to the wall for the neutral species, (r„ ■ n)(.^ = 0 , L) = |y/8T,i/7rm„ with n 
the normal vector to the corresponding electrode, and from the assumption of the electric 
field dominated particle flux to the wall for the ions, (F, ■ n)( 2 ; = 0, L) = a/XjnjE, with a the 
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switching function, taking value 1 when E ■ n > 0 and 0 when E ■ n < 0, respectively (e.g. 

[I])- 

It is essential to correctly determine the momentum exchange frequency i'hs,N in Eq- ([^ 
and the mobilities in Eq. ([^ (note that one of these quantities can be calculated from the 
other utilizing Eq. (|^) caused by collisions of ion species with the helium background gas. 
Although ion mobilities are often assumed to be constant (see, e.g., ID), a better approxi¬ 
mation attempted from the physics considerations yielding dependence of the mobilities on 
the electric field demonstrates that such a dependence can matter [TH] . Following the latter 
argument, in the present work we propose to calculate mobilities and the corresponding 
momentum exchange frequency with help of the following auxiliary Monte-Carlo code. 

The code resembles very much a PIC/MCC code where only one of the ion species is 
traced. The particles are evolved in time under action of a prescribed constant electric held 
Eq and collisions with the background neutral gas having a temperature Tn. The collisions 
are modeled using the Monte-Carlo method with the same energy resolved cross-sections 
as in the full PIC code. At the initial moment the particle velocities are sampled from 
a maxwellian distribution with an initial temperature T^o. The computational domain is 
assumed to have periodic boundaries. After several nanoseconds the ions acquire the sta¬ 
tionary drift velocity (see Figure 6 in HT] and the accompanying discussion in that reference), 
which can be determined by calculating the total mean velocity of the particle ensemble. 
Carrying out this procedure for a number of the electric held values, one can construct a 
look-up table or an analytic ht of mobility (calculated as u/Eq, where u is the drift velocity) 
and momentum exchange frequency (calculated from Eq. ([^) versus the electric held and 
the drift velocity, respectively. Following the suggestion in w, we use both in the full 
PIC and the Monte-Carlo code for the ion isotropic scattering the cross-section equal to 
Qm = Qi + ‘^Qb, where Qi is the elastic isotropic cross-section and Qh is the elastic backward 
scattering cross-section. 

Note that this technique is more general than the one suggested in [TH]. It enables an 
accurate calculation of the momentum exchange frequencies and mobilities for any species, 
for which energy-resolved collision cross-sections are known, and under any conditions of 
practical interest (for example, when the colliding species have diherent velocity distribution 
functions with comparable but distinct characteristic energy/temperature). Calculation of 
the corresponding collision frequencies and mobilities has to be performed only once for a 


10 


given pressure of the background helium gas and thus does not cause any computational 
overhead in the hybrid code simulations. 

D. Coupling Between the Electron and Ion Models 

Coupling of the kinetic electron and the fluid ion models occurs through the net charge 
density in the Poisson’s equation, which includes both electron and ion densities and through 
any reaction involving both electrons and ions. Among the typical reactions are the sec¬ 
ondary electron emission caused by ions impinging on the electrode surfaces, production 
of electrons and positive ions in the electron impact. Penning and the metastable pooling 
ionization reactions, and recombination reactions of electrons with positive ions. 

Once the charge densities are calculated, the Poisson’s equation has to be solved. In order 
to achieve a needed regime of operation for the discharge in numerical simulations one has 
to limit the current flowing through the discharge. If it is not done (for example, if a hxed 
voltage amplitude RF source is used in the simulations), the number of superparticles and 
hence the particle density tend to diverge with time as the modeled discharge spontaneously 
goes into the 7 regime with an ever increasing current. In order to limit the total current 
in the simulations one can chose either to use a hxed amplitude currrent source or, if a 
hxed amplitude voltage source is used, either to use an external resistance or to limit the 
power absorbed in the discharge by adjusting the voltage amplitude to meet the prescribed 
power. The latter approach seems to correspond to experimental observations better, the 
latter showing generation of additional harmonics rather in the measured current than in 
the voltage signal |1]. To adapt the voltage ampltide for matching the prescribed power 
we have implemented the following simple algorithm. Starting from an initial guess for the 
voltage amplitude [/°, during each RF period the voltage ampltiude remains hxed and the 
period-averaged power absorbed by the discharge during the A^th period is calculated as 
^ V{j + eodE/dt)dt, where S is the electrode area, V is the voltage, T is the 

RF period, j and eodE/dt are the conduction current and the displacement current densities 
at the driven electrode. The new value for the voltage amplitude during the {N + l)th RF 
period is calculated as 

pN+l ( 11 ) 

with Pq the prescribed power and a some numerical parameter, which can be chosen to be 
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small for a smoother or large for a faster adjustment, respectively. It can be also be chosen 
to be time-adaptive, for example, being larger at the beginning of the simulation and smaller 
as the simulation proceeds. 

The secondary electron emission algorithm is implemented by employing the following 
technique. The number of heavy particles (ions or metastables) ANi hitting an electrode 
during a time interval At equals to TiSAt. Each such heavy particle can produce a sec¬ 
ondary electron with probabilty equal to 7 j, which results in an average number of sec¬ 
ondary electrons produced during a time step equal to Nsec.ei = liANi. This translates 
to Nsec.ei/Afp/sp superparticles with Np/gp dehned in Eq. ([^. Thus, during each time step 
one produces int {Nsec.ei/Afp/sp) (with int representing the integral part of a number) elec¬ 
trons and performs a comparison of the R uniformly distributed pseudorandom number with 
Nsec.ei/Np/sp — int [Nsec.ei/Np/sp). If the former is smaller than the latter, one additional sec¬ 
ondary electron superparticle is produced. Such a Monte-Carlo technique ensures correct 
number of secondary electrons generated on average. 

The electron impact ionization is treated according to the modihed null-method described 
in m When an impact ionization or a metastable excitation event occurs, a corresponding 
location of the newly produced ions or metastable species are recorded and the corresponding 
ion or metastable density is incremented by Np/gp/AV with AV the grid cell volume along 
with creation of a new electron superparticle. Velocity of the ejected electrode is determined 
by a random scattering on a sphere in velocity space with radius corresponding to the energy 
of the ejected electron [T5] . 

Ionization events due to the Penning or the metastable pooling reactions are treated 
in a similar fashion, but with the number of electron superparticles corrected with regard 
to the corresponding reaction. For example, for the Penning ionization such number is 
equal to ANe/Np/gp, where number of electrons produced during a time step is equal to 
ANe = Rn^.n]siAtAzS with R the corresponding reaction rate, n* the density of the cor¬ 
responding metastable species and njsi the density of the background neutral gas ionized 
by the metastable species. As in the case of the impact ionization, velocity of the electron 
created in such a process is to be found by scattering on a sphere in velocity space with 
radius calculated from the corresponding energy (different for each of the processes). 

The recombination reactions between electrons and positive ions are accounted for numer¬ 
ically as follows. First, the maximum collision probability during a time step is calculated 
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TABLE I: Electron-neutral collisions 


(Rl) e -I- He — )• e -|- He ISO] 

(R2) e-l-He—>-e-|-e-|- He"'' (t(T) [20] 

(R3) e + He ^ e + He* a{8) [50] 

(R4) He^ -|- He — )■ He"'' + He Qi{£) + 2Qb{£) [22], see also text 

as Pmax = 1 — exp(—with Umax = max?7,j(^) ■ max(cr(£^)v). Then, the rest of the 

2 £ 

null collision procedure goes as the previously mentioned modihed null-collision method of 
[HI with a correction respecting the nonuniform ion density profile when calculating the 
actual collision probability. When a recombination event occurs, the corresponding electron 
superparticle is removed from the computational domain and the corresponding ion density 
is decremented by Np/sp/AV. 

III. VALIDATION OF THE HYBRID CODE 

The approach described in the previous section enables an efficient implementation of 
numerical model for discharges posessing a complicated chemistry. The goal of the present 
work is, however, to validate the approach itself and to demonstrate its applicability to 
adequate modeling of the highly collisional plasma discharges. To this end in the next 
section we compare the simulation results of the hybrid code with the simulation results of 
the GPU PIC/MCC code. The latter has been verified in the benchmark study of [16], but 
this time it has been ran assuming the atmospheric pressure in a pure He discharge using 
a simple reduced chemistry set similar to the one used in [S] and listed in Table Such a 
chemistry set does not represent the actual physics taking place in such discharges because 
metastable atoms and excimer molecules are neglected, whereas they usually dominate the 
electron production through the pooling reactions. However, the simplihed chemistry allows 
to focus on validating if the hybrid approach properly describes the particle, momentum, and 
energy transport, which should match those of the PIC code making minimum assumptions 
and thus describing the physics in such discharges comprehensively. Hence, in this reduced 
set only electrons and He+ ions are tracked. 

For a pure helium discharge we took the energy resolved elastic cross-sections for the 
collisions between the electrons and the He background neutral gas from the XPDPl code 
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TABLE II: Mobilities, [mV(Vs)] 

1.18 - 0.185(|E|/E;o) , \E\ < Eq = 1.47 x 10® V/m 
0.94 (Eo/|E|)0'3 , \E\ > Eo 


TABLE III: Momentum Exchange Erequency, 10^® [s 
2.07 + 0.88(|'u|/«o)^-®® , 1^1 < uo = 2 X 10® m/s 
1.54 + 1.36(|u|/tto) , |tt| > uq 


[ 20 ] (for the benchmarking pnrposes a particular choice of these cross-sections does not 
play a signihcant role, since the electron component is modeled using the same approach 
and the same parameters, including the cross-sections, both in the hybrid code and in 
the fully kinetic GPU PIC/MCC code), and for the collisions of He"*" ions with the He 
neutrals from [ 22 ], Qi{S) = 7.63 x [m^] and Qb{S) = 10 “^®/[(£ 1 / 1000)®-^®(1 -|- 

£^/1000)®'^®(l -|- [m^] with S the relative energy in eV calculated in the center of 

mass Mrvl^i/2 ^ MneVrei/^^ where Mr is the reduced mass of the collision partners and 
Vrei is their relative velocity. Utilizing these cross-sections, the Monte-Carlo code described 
in the previous section yielded mobilities for the He"*" ions in a He gas under atmospheric 
pressure, which we approximated by the analytical hts listed in the Table A better £t 
can be made to match expected the high held asymptotics, fi oc U”® ® (see the discussion in 
[T 8 ] ), more accurately, yet we have found that the simple formula given in Table gives a 
good agreement for E < 5x 10® V/m, which is the range of the electric held values observed 
in the simulations. Similarly, the analytical £t for the momentum exchange frequency in 


the elastic scattering between He"*" ions and He atoms is listed in Table The secondary 
electron emission coefficients were calculated using the empirical formula given in [ 23 ] (the 
reference quotes 50% accuracy for this formula), ~ 0.016{Eiz—2E^) with Ei^ the ionization 
potential of the incident ion and E^ ^ 4.5 eV the work function. This yields 7 ~ 0.25 for 
the helium ions. 

For the hrst set of benchmarks testing the hybrid code against the particle-in-cell code 
for the atmospheric pressure capacitively coupled plasma discharge we have chosen to study 
two different regimes of the discharge operation, an ohmically heated discharge (H mode 
with Pq/S = 10"^ W/m^ and I = 1.5 mm (Case I) and a discharge dominated by the 
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FIG. 1: Shown are RF period-averaged density profiles for the Case I (11 mode) calculated with 
the PIC code and the different versions of the hybrid code, see text. 

ionizing avalanches produced by the secondary electrons accelerated by strong electric held 
in a discharge with the electrode distance I = 75 /rm comparable to its sheath width, so 
that we call it a ’’short” 7 discharge (Case II). For the latter discharge the power density is 
considerably higher, we studied a case with Po/S = 2.4 x 10® W/m^. 

Fig. 1 shows simulation results of the Case I for the RF period-averaged particle denisty 
prohles provided by the PIC/MCC code and three different versions of the hybrid code. The 
hrst version of the hybrid code (”HC1”) employs the full fluid model described in the previous 
section, which accounts for the ion inertia by solving Eq. (|^ and uses the analytic £t for the 
momentum exchange frequency given in Talbe |T], determined with help of the Monte-Carlo 
code described in the previous section. The second version of the hybrid code (” HC2”) uses 
the drift-diffusion approximation with the analytic £t for the He’*' mobilities in He given 
in Table also obtained from the results of the auxiliary Monte-Carlo code. Finally, the 
third hybrid code version (”HC3”) uses the drift-diffusion approximation described in the 
previous section and a simple constant approximation for the mobility value (taken from 
the same analytic fit at zero electric field, which incidentally turns out to be very close to 
the value adopted in |1]). The comparison between simulation data reveals that this regime 
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FIG. 2: The RF period-averaged density profiles for the Case II (7 mode, ions with Xe mass) 
obtained by the PIC/MCC code and different versions of the hybrid code, see text. 

is quite accurately described by all the considered fluid models. As expected, the largest 
deviation from the kinetic results is demonstrated by the constant-mobility model (HC3), 
yet its results are still very close to prediction of the PIC/MCC code. 

Fig. 2 shows results of the benchmark of the PIC/MCC code and the same versions of 
the hybrid code as described above used for parameters of the Case II (7 mode, in which 
the discharge is sustained by the ion-induced secondary electron emission and the ionization 
avalanches caused by the secondary electrons accelerated by the electric held). In this regime 
plasma density is much higher than in the G mode and the bulk plasma is quasineutral on 
average. Therefore, we have plotted only the RF period-averaged electron densities and 
omitted the ion densities, which are equal to the electron densities in the major part of the 
discharge, in order to simplify the hgure. One can see that the drift-diffusion approximation 
gives a relatively good agreement with the prediction of the PIC code and the full huid model 
used in the HCl version of the hybrid code. The three prohles (PIC, HCl, and HC2) are 
very close over almost entire gap and differ only slightly in the center region. In contrast, the 
hybrid code version using the constant mobility (the zero electric held value) demonstrates 
an average deviation from the kinetic result. The difference between the kinetic result and 
the result of the constant mobility drift-diffusion model can be attributed to the fact that 
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TABLE IV: Mobilities, [mV(Vs)] 
3.05 - 0.75(|E|/Eo) , \E\ < Eq = 2 x 10® V/m 
2.3 (Eo/|E|)0-37 , \E\ > Eo 


TABLE V: Momentum Exchange Frequency, 10® 
6 + 7 X 10“®|up'®^ , |u| < Mo = 4 X 10^ m/s 
4.5 + 27|m|/mo , |m| > mq 


electric field in the studied example is much higher than in the 12 mode and thus dependence 
of the helium ion mobility on the electric field given in Table [IT] starts to play a role. Still, 
the discrepancy between the peak values obtained with the PIC and the HC3 codes are less 
than 4% and thus are not significant. 

In [T7j it was suggested that heavier ions should require more time to reach the stationary 
drift velocity under influence of electric field and collisions with the background gas. To see 
how ion mass affects the accuracy of the fluid models used in the hybrid code we have 
performed the following test. We have used the same ion-neutral collision crosssections as 
for the He^ - He collisions given in the Table before, but this time we have increased the 
ion mass to match that of xenon. The Monte-Carlo code used to calculate the momentum 
exchange frequency and mobilities has provided data, which we fitted analytically as shown 


in Tables IV and |Vj The corresponding simulations conducted with parameters similar to 
the Case I (H mode) have shown very small discrepancy between the models and is not 
shown here. However, an analog of the Case H (’’short” 7 discharge) simulated with the 
heavier ions (henceforth called Case HI) indeed exhibits dependence of the resulting RF- 
averaged electron density profiles on a particular model (see Fig. 3). Whereas the profiles 
calculated with the PIC/MCC code and the full fluid model (HCl) are still quite close to 
each other, the drift-diffusion models with variable mobility (HC2) and constant mobility 
(HC3) demonstrate a 20% digression from the kinetic result. Interestingly, the latter models 
yield different signs of the deviation from the kinetic result. 

Finally, to test how accurate are the ion fluid models at lower pressures we have studied 
an example of a pure He discharge operated at 300 mTorr pressure of the He working gas, 
Pq/S = 170 W/m^ and I = 3 cm. The corresponding momentum exchange frequencies and 
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FIG. 3: The RF period-averaged density profiles for the Case III (7 mode) obtained by the 
PIC/MCC code and different versions of the hybrid code, see text. 


TABLE VI: Mobilities, [mV(Vs)] 

2.2 - 0.78(|E|/Eo) , \E\ < Eq = 2.8 x 10^ V/m 
1.42 (Eo/E)0'41 , \E\ > Eo 


mobilities to be used in the ion fluid models were calculated with the Monte-Carlo code 
described before and £t with analytic functions as shown in Tables VI and VII[ At such a 
pressure the momentum exchange frequency becomes comparable to the driving frequency 
and one can expect that the explicit ion flux time modulation on the left hand side in 
Eq. Q becomes comparable with the right hand side. This should lead to a breakdown of 
the drift-diffusion approximation. Indeed, the corresponding results (see Fig. 4) demonstrate 
that even the G mode (more commonly referred to as the a mode at such pressure) is not 
properly simulated by the hybrid code using the drift-diffusion simulation with a realistic 
mobility (HC2). In contrast, the hybrid code version utilizing the full ion fluid model is in 
an excellent agreement with the kinetic result. Therefore, despite the drift-diffusion model is 


TABLE VIL Momentum Exchange Frequency, 10® [s 
8 -|- 0.49|M|/2.e4 , u < 10® m/s 
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FIG. 4; The RF period-averaged density profiles for the Case IV (11(a) mode, 300 mTorr) obtained 
by the PIC/MCC code and different versions of the hybrid code, see text. 

frequently used in the literature also at low pressures, one must be very careful in interpreting 
its results, as its assumptions are very likely to be violated there. 


IV. CONCLUSIONS 

The present work describes a hybrid numerical scheme that can be used for simulations of 
highly collisional discharges with a complex chemistry. The scheme uses a kinetic description 
for electrons based on the PIC/MCC method and considers several possible fluid models 
for description of ion species. For the ’’full” fluid model accounting for the explicit time 
modulation of the ion flux a simple numerical scheme is proposed. Its results are confronted 
with results of the purely kinetic PIC/MCC code and the popular drift-diffusion approach 
for several exemplary discharges in regimes of practical interest. It is demonstrated that the 
drift-diffusion model with the constant mobility performs well for plasma discharges under 
atmospheric pressure unless the discharge is operated deeply in the 7 regime with heavy ion 
species present. It is also shown that the drift-diffusion approximation breaks down at lower 
pressures, whereas the hybrid code version with the ’’full” fluid model remains very close to 
predictions of the kinetic code. 

It is worth noting that the proposed ion fluid model taking into account ion inertia imposes 
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only a minor complication of the numerical algorithm compared to the popular numerical 
schemes using the drift-diffusion approximation, which consists in solving an additional 
equation for the ion flux. However, considering the substantial improvement in the accuracy 
of the physics description of the former, we suggest that it should generally be preferred 
against the latter, at least as far as RF-driven discharges are concerned. 
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